14th  International  Conference  on  Information  Fusion 
Chicago,  Illinois,  USA,  July  5-8,  2011 


Forecasting  Probability  of  Target  Presence  for 
Ping  Control  in  Multistatic  Sonar  Networks 
using  Detection  and  Tracking  Models 


Cherry  Y.  Wakayama 

Code  56560 

SPAWAR  Systems  Center  Pacific 
San  Diego,  CA,  U.S.A. 
cherry.wakayama@navy.mil 


Doug  J.  Grimmett 

Code  56560 

SPAWAR  Systems  Center  Pacific 
San  Diego,  CA,  U.S.A. 
doug .  grimmett  @  navy.mil 


Zelda  B.  Zabinsky 

Industrial  &  Systems  Engineering 
University  of  Washington 
Seattle,  WA,  U.S.A. 
zelda@u.  washington.edu 


Abstract — This  paper  describes  the  forecasting  of  probability 
of  target  presence  in  a  search  area  (also  referred  to  as  the  PT 
map)  considering  both  detection  and  non-detection  conditions. 
Tracking  results  are  also  incorporated  to  obtain  a  more  accurate 
PT  map  under  the  detection  condition.  The  probability  of  target 
presence  is  a  suitable  metric  for  real-time  ping  control  for  a 
submarine  search  mission,  whose  objective  is  to  quickly  identify 
and  localize  as  many  targets  as  possible  within  the  search 
area.  Existing  formulations  of  the  probability  of  target  presence 
metric  for  ping  control  include  an  open-loop  approach  in  which 
measurements  are  ignored  or  a  semi-adaptive  approach  in  which 
measurements  are  considered  but  without  the  true/false  target 
investigation.  Since  false  contacts  are  inevitable  in  practical 
applications  and  the  true/false  target  investigation  of  the  contacts 
is  not  immediate,  tracking  results  must  be  considered  in  the  PT 
map  generation  to  obtain  an  accurate  assessment  of  the  present 
and  projected  operational  pictures.  We  develop  an  approach  to 
obtain  the  current  and  forecasted  PT  maps  by  incorporating  a 
measurement  model,  a  sonar  performance  model,  Bayes  theorem 
and  a  centralized  Kalman-Filter  based  tracker.  The  PT  map  is 
composed  of  two  portions:  the  portion  which  contains  detected 
target  probability  and  the  portion  which  contains  missed  target 
probability.  Each  portion  of  the  PT  map  is  updated  and  propa¬ 
gated  separately.  The  forecasted  PT  map  at  the  next  ping  time 
is  obtained  by  combining  the  two  propagated  PT  maps.  It  will 
be  demonstrated  by  simulations  that  the  combined  forecasted  PT 
map  represents  an  accurate  multistatic  operational  picture  and 
can  be  used  with  a  sonar  performance  model  to  obtain  a  field 
metric  for  ping  control  optimization  for  the  area  search  mission. 
Keywords:  Ping  control,  sensor  management,  multistatic 
tracking. 

I.  Introduction 

Multistatic  active  sonar  networks  provide  an  Anti- 
Submarine  Warfare  (ASW)  capability  against  small,  quiet, 
threat  submarines  in  the  harsh  clutter- saturated  littoral  and 
deeper  ocean  environments  through  the  expanded  geometric 
diversity  of  a  distributed  field  of  sources  and  receivers.  How¬ 
ever,  such  networks  cannot  automatically  exploit  their  full 
potential  without  intelligent  management  and  control  methods 
for  the  sources  and  receivers  because  of  the  variability  in 
acoustic  environmental  conditions,  sensor  performance  and 
threat  target  behavior  [11.  In  general,  control  methods  and 


tactical  decision  aides  may  be  applied  to  address  the  sen¬ 
sor  placement,  signal  and  information  processing  or  sonar 
ping  control  problems.  In  addition,  the  control  methods  and 
metrics  for  dealing  with  threat  submarines  may  be  mission 
dependent  in  order  to  effectively  optimize  the  sensor  network 
performance  for  the  given  mission  objective(s).  This  paper  is 
concerned  with  the  sensor  ping  control  method  for  the  area 
search  mission  type,  whose  objective  is  to  quickly  identify 
and  localize  as  many  threat  targets  as  possible.  In  this  case, 
the  probability  of  target  presence  is  an  effective  metric,  and 
the  adaptive  control  approach,  in  which  the  evaluation  and 
optimization  of  the  metric  make  use  of  all  the  information 
obtained  during  the  operation  to  enhance  real-time  target  De¬ 
tection,  Classification,  and  Localization  (DCL)  performance,  is 
appropriate.  In  other  cases,  the  mission  may  be  to  ensure  that 
some  selected  area  is  free  of  threat  targets.  As  discussed  in  [2], 
the  ASW  Residual  Risk  metric  is  an  effective  metric  for  the 
area  clearance  mission  type  and  the  sensor  control  approach 
may  be  optimized  prior  to  the  operation. 

Approaches  to  evaluating  the  probability  of  target  presence 
metric  for  ping  control  solutions  for  the  area  search  mission 
are  presented  in  [3]  and  [4].  The  formulation  presented  in 
[3]  ignore  the  measurements  and  assumes  that  no  detection 
has  been  confirmed.  The  probability  of  target  presence  in 
each  xy  grid  cell  within  the  search  area  (also  referred  to 
as  the  PT  map)  is  updated  in  a  Bayesian  manner  under 
the  non-detection  hypothesis.  In  [4],  the  measurements  are 
considered  in  the  formulation  and  the  PT  map  is  updated  in 
a  Bayesian  manner  under  the  detection  hypothesis  if  there  is 
a  measurement  in  the  xy  grid  cell  or  under  the  non-detection 
hypothesis  otherwise.  In  both  formulations,  the  forecasted  PT 
map  is  obtained  by  propagating  the  posterior  probability  using 
a  diffusion  process  model.  In  the  presence  of  false  contacts, 
which  is  inevitable  in  practical  applications,  one  must  also 
account  for  the  track  consistency  in  order  to  discriminate  the 
target-originated  contacts  against  the  contacts  that  are  likely  to 
be  false  alarms  for  obtaining  accurate  multistatic  operational 
pictures.  The  true/false  target  investigation  is  not  considered 
in  the  existing  formulations. 
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We  develop  an  approach  to  construct  the  PT  maps  which 
considers  the  measurements  (sonar  contacts)  as  well  as  the 
true/false  target  investigation  of  the  contacts  (via  tracking)  and 
evaluates  true  and  false  target  probabilities,  missed  target  prob¬ 
abilities  and  no  target  probabilities  by  incorporating  a  sonar 
performance  model,  a  measurement  model,  Bayes  theorem, 
and  a  multistatic  target  tracker  (MTT).  We  assume  that  an 
xy-grid-based  probability  of  detection  map  for  each  source- 
receiver-waveform  combination  is  provided  by  a  sonar  per¬ 
formance  model.  Our  objective  is  to  obtain  xy-grid-based  PT 
maps  at  the  current  and  next  ping  times  from  the  measurements 
and  the  tracker  outputs.  In  addition  to  the  PT  maps,  it  is  also 
required  to  obtain  the  multiple-target  state  estimates.  While 
performing  the  area  search  mission,  it  is  desired  to  track  not 
only  the  possible  target-originated  contacts  but  also  the  area 
that  has  been  not  been  searched.  The  complete  PT  maps  should 
contain  information  resulting  from  the  presence  and  absence 
of  contacts  in  the  area  that  has  been  searched. 

A  portion  of  the  probability  of  target  presence  in  each  grid 
cell  due  to  the  possible  target-originated  contacts  is  obtained 
from  the  multiple-target  state  estimates  of  the  tracker.  The 
missed  detection  probability  in  the  area  that  has  been  searched 
is  obtained  by  incorporating  a  sonar  performance  model,  a 
measurement  model  and  Bayes  theorem.  The  tracker  in  this 
work  is  a  conventional  centralized  Kalman-Filter  (KF)  based 
tracker  which  performs  explicit  data  associations  between 
measurements  and  targets  to  update  the  multiple-target  state 
estimates  [5].  This  tracker  does  not  assume  any  a  priori 
distribution  of  targets  and  does  not  account  for  the  missed 
targets  until  tentative  tracks  have  been  initiated. 

Alternatively,  the  PT  maps  may  also  be  generated  from 
the  implementations  of  Probability  Hypothesis  Density  (PHD) 
filter  [6]  and  [7].  The  PHD  filter  operates  on  the  single-target 
state  space  and  avoids  the  explicit  data  association  between 
measurements  and  targets.  It  is  a  recursion  that  propagates 
the  first-order  statistical  moment,  or  intensity,  of  the  random 
finite  set  of  states  according  to  given  model  parameters  that 
define  target  survival  probability  and  the  shapes  of  target 
birth  and  spawning  intensities  [7]  and  [8].  In  order  to  obtain 
the  multiple-target  state  estimates,  one  must  then  perform 
peak  extraction  from  the  posterior  intensity.  In  this  paper,  we 
assume  the  use  of  a  conventional  KF  based  tracker  and  develop 
a  formulation  for  the  PT  map. 

This  paper  proceeds  as  follows.  In  Section  II,  we  describe 
the  formulation  of  the  probability  of  target  presence  forecaster. 
We  present  simulation  results  to  demonstrate  the  performance 
of  the  forecaster  in  Section  III.  Conclusions  and  future  work 
are  provided  in  Section  IV. 

II.  Formulation  of  Probability  of  Target 
Presence  Forecaster 

The  objective  of  the  forecaster  is  to  obtain  accurate  oper¬ 
ational  pictures  of  the  multistatic  network  for  real-time  ping 
control  by  evaluating  true  and  false  target  probabilities,  missed 
target  probabilities  as  well  as  no  target  probabilities.  This  is 
achieved  by  updating  PT  maps  after  each  ping  using  Bayes 


theorem  under  the  detection  and  non-detection  hypotheses  sep¬ 
arately.  The  tracking  (true/false  target  investigation)  results  are 
incorporated  to  further  refine  the  PT  map  under  the  detection 
hypothesis.  We  refer  to  the  two  PT  maps  generated  under 
the  detection  and  non-detection  hypotheses  as  the  positive 
PT  map  and  the  negative  PT  map,  respectively,  for  the  ease 
of  discussion.  The  positive  PT  map  describes  the  detected 
target  probability  and  the  negative  PT  map  describes  the 
missed  target  probability.  The  mathematical  description  of  the 
forecaster  is  given  in  the  following. 

We  define  a  search  region  on  a  xy  coordinate  space.  We 
discretize  the  search  region  into  a  grid  G  =  [gi,g2,mmm  ?  9  n9]t 
where  =  {{x,y),x  e  [x*  -  <5/2,  x*  +  <5/2 ),y  e  [yt  - 
S/2,  yi  +  <5/2)},  Ng  is  the  total  number  of  grid  cells,  (x*,  y, ) 
represents  the  center  coordinates  of  cell  i,  and  S  is  the  width 
of  each  cell.  The  state  vector  at  time  tk  on  G  is  represented  by 
T{k)  =  [T1(k),T2(k),  •  •  •  ,TNg(k)]T ,  where  Tt(k)  is  a  binary 
random  variable:  Ti(k)  —  1  represents  the  event  that  at  least 
one  target  is  in  cell  i  at  time  tk,  and  Ti(k)  =  0  represents  the 
event  that  no  target  is  in  cell  i.  The  PT  map  at  time  tk  is  de¬ 
fined  as  [P(Ti(*0  =  l),P(T2(fc)  =  1),--;  ,P (TNg(k)  =  1)]T 
Note  that  the  summation  of  the  probabilities  on  the  PT  map 
represents  the  expected  number  of  targets  in  the  search  region. 

We  denote  the  set  of  measurement  contacts  at  time  tk  by 
Z{k)  —  { Zx(k ),  Z2(k),-  ■  •  ,  ZnR(k)},  where  tir  is  the  num¬ 
ber  of  receivers.  The  measurement  set  at  receiver  j  at  time  tk  is 
represented  by  Z^(k)  —  { Z{(k ),  ZJ2(k ),  •  •  •  ,  Z3n.  (k)},  where 
rij  is  the  total  number  of  contacts  at  receiver  j.  We  model 
each  measurement  variable  Z Jq  as  a  Gaussian  distribution  with 
a  mean  vector  of  [xJq,yJq]T  and  a  covariance  matrix  of  RJq 
for  q  =  1,  2,  •  •  •  ,  rij  and  j  =  1,  2,  •  •  •  ,  tir.  The  mean  vector 
and  covariance  matrix  are  derived  from  the  time-of-arrival  and 
bearing  measurements  of  contact  q  at  receiver  j  and  their 
associated  measurement  errors  as  described  in  [9]. 

We  define  an  observation  vector  at  time  tk  on  G ,  which  is 
denoted  by  y(k)  =  [Yi(/c),  Y2{k)f  •  ,  YNg  (k)]T,  where  Y^k) 
is  a  binary  random  variable:  Yi(k)  —  1  represents  the  event 
that  measurement  observation  occurs  in  cell  i ,  and  Yi(k)  =  0 
represents  the  event  that  no  measurement  observation  is  in  cell 
i.  We  refer  to  the  event  Yi  =  1  as  the  detection  hypothesis  and 
Yi  =  0  as  the  non-detection  hypothesis. 

A.  Bayes  Update 

For  the  ping  generated  at  time  tk,  we  receive  rij  contacts 
at  receiver  j.  We  are  going  to  suppress  k  for  ease  of  notation. 
The  probability  that  the  measurement  data  Z Jq  occurs  in  cell  i 
is  given  by: 

P((1Z)^  =  l\Z{)  =  [  N{[x{,yi\T,R{)dxdy,  (1) 

where  Af(m,  E)  denotes  a  Gaussian  density  with  mean  m 
and  covariance  E,  and  (Yi)Jq  =  1  represents  the  detection 
hypothesis  with  respect  to  measurement  q  of  receiver  j. 
Aggregating  the  probability  information  across  all  contacts 
at  receiver  j ,  the  probability  that  there  is  a  measurement 
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observation  in  cell  i  is  obtained  as: 


P(1^J  =  1|ZJ)  —  min 


i.Ep((y*)«  =  iiz«) 


(2) 


\  q=  1  / 

The  probability  that  there  is  no  measurement  observation  in 
cell  i  after  considering  the  measurement  data  ZJ  at  receiver 
j  is  given  by: 


P (Y?  =  0| Z*)  =  1  -  P (Yi  =  l\Zi).  (3) 


The  grid  size  should  be  chosen  sensibly  with  respect  to 
measurement  density  such  that  the  summation  of  probabilities 
in  Eq.  (2)  will  not  exceed  1  in  most  cases. 

The  probabilities  of  target  presence  in  each  grid  cell  un¬ 
der  detection  and  non-detection  hypotheses  at  receiver  j  are 
obtained  using  Bayes  theorem  as: 


P(T,  =  1|17  =  1) 

PCTi  =  i|y?  =  o) 


PD]  •  PT* 

PD]  •  PT*  +  Pfa]  (1  -  PT*)  ’ 

_ (1  -  PD])PT* 

(1  —  PD])PT*  +  (1  —  Pfa])(l 


(4) 


PT  i)' 

(5) 


where  PT*  =  P(T*  =  1)  represents  the  a  priori  probability  of 
target  in  the  grid  cell  i,  PD]  =  P(Y^  =  1|T*  =  1)  represents 
the  probability  of  detection  at  receiver  j  given  there  is  a  target 
in  the  grid  cell  i,  and  Pfa]  =  P(Y^'7  =  1|T*  =  0)  represents 
the  probability  of  false  alarm  at  receiver  j  in  the  grid  cell  i. 
The  PD  map  for  each  source-receiver-waveform  combination 
is  precomputed  using  a  sonar  performance  model  which  will 
be  discussed  later. 

The  PT  map  for  receiver  j  after  updating  the  measurement 
data  Zi  is  obtained  as: 


P(T*  =  1|  Z*)  =P  (T*  =  1\YP  =  1)P  (YP  =  1|  ZP) 

+  P(T*  =  1|  YP  =  0)P  (YP  =  0|  Zj),  (6) 
=P(T*  =  l\ZP)*  +  P(Ti  =  l\Zi)~y  (7) 

where  P (T*  =  1| ZJ)+  and  P(T*  =  1| Z^)~  represents  the 
positive  and  negative  PT  maps  for  receiver  j  respectively. 
The  Ur  positive  PT  maps  generated  by  tir  receivers  are 
then  fused  to  obtain  a  single  positive  PT  map  after  each 
ping;  a  negative  PT  map  is  obtained  similarly.  Some  source- 
receiver-waveform  combinations  provide  more  confidence  in 
their  detections  at  certain  grid  cells  according  to  their  PD 
values  given  by  sonar  performance  model  predictions.  We 
define  a  weighting  coefficient  for  receiver  j  at  grid  cell  i  as: 


PD] 


YT^pdI 


(8) 


The  fused  positive  and  negative  PT  maps  are  then  obtained  as 
the  weighted  average  of  the  respective  PT  maps  given  by: 


nR 

P(Ti  =  l|Z)+=^W|P(Ti  =  l|^)+,  (9) 

3  = 1 

riR 

p  (Ti  =  1|  Z)~  =J2wip  (Ti  =  (10) 

3  = 1 


We  update  the  PT  maps  separately  under  the  detection 
and  non-detection  hypotheses,  because  we  further  refine  the 
positive  PT  map  using  the  tracking  results  which  will  be 
discussed  in  the  next  subsection.  We  propagate  the  two  PT 
maps  separately  and  combine  them  after  propagation  to  obtain 
the  forecasted  PT  map. 

The  values  of  PD]  and  Pfa]  are  obtained  from  a  sonar 
performance  model  for  each  source-receiver-waveform  com¬ 
bination  with  an  average  omni-directional  target  strength  in 
the  search  area  [10].  Assuming  the  layout  of  the  sources 
and  receivers  is  fixed,  the  average  PD  maps  for  different 
source-receiver-waveform  combinations  can  be  precomputed 
and  stored.  We  assume  that  the  average  false  alarm  rate  is 
constant  and  known  for  each  receiver;  Pfa]  is  computed  from 
this  given  false  alarm  rate,  the  area  of  the  search  region,  and 
the  area  of  the  grid  cell.  We  also  assume  that  at  the  beginning 
of  the  scenario  each  grid  cell  has  a  target  probability  of  p ,  i.e. 
P{Ti( 0)  =  1)  =  p.  After  each  ping,  we  generate  a  forecasted 
PT  map  for  the  next  ping  time,  which  becomes  the  a  priori 
PT  map  for  processing  the  next  ping  measurement  data. 

B.  Multi-Target  Tracker 

The  tracker  processes  contacts  and  identifies  possible  targets 
in  the  presence  of  false  alarms.  This  task  requires  the  removal 
of  a  large  number  of  false  contacts,  the  association  of  true 
contacts,  and  the  fusion  of  contact  information  to  estimate 
the  location  and  velocity  of  targets  of  interest.  There  is  much 
literature  on  sensor  data  fusion  and  target  tracking;  classical 
references  include  [11]  and  [12].  Various  types  of  trackers  can 
be  implemented.  However  an  appropriate  tracker  implementa¬ 
tion  must  consider  the  tradeoff  between  complexity  and  quality 
depending  on  the  operation  scenarios.  In  this  work,  we  use  a 
centralized,  Kalman  Filter  (KF)  tracker  described  in  [5]  for 
the  real-time  adaptive  ping  control  application.  Although  the 
current  MTT  implementation  uses  a  naive  nearest  neighbor 
data  association  scheme,  it  will  serve  as  a  baseline  tracker 
to  illustrate  the  procedure  of  incorporating  the  KF  estimates 
to  the  PT  maps.  The  performance  of  the  forecaster  can 
certainly  be  improved  by  implementing  a  more  sophisticated 
data  association  scheme  such  as  Multi-Hypothesis  Tracking 
(MHT)  at  the  cost  of  computational  complexity  [12]. 

We  briefly  describe  the  tracker  implementation  in  the 
following.  The  input  to  the  tracking  algorithm  is  a  series  of 
contact  files  (scans),  unique  to  each  source-receiver- waveform 
and  time  of  ping  transmission,  provided  by  the  sensor 
network.  In  the  current  formulation,  we  consider  only  FM 
waveforms  with  a  constant  ping  interval. 

Target  Motion  Model 

The  target  state  vector,  X ,  consists  of  the  x  and  y 
components  of  both  position  and  velocity,  i.e.  X(t)  = 
[x(t),  y(t),  x(t),  y(t)]T.  The  discrete-time  nearly  constant  ve¬ 
locity  motion  model  is  given  by: 


^c+i  =  &kXk  +  cjfc,  Wk  ~  A/*(0,  Qk)  (11) 
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where  X&+i  is  the  target  state  at  time  tk+i,  &k  is  the 
state  transition  matrix,  and  Qk  is  the  covariance  matrix 
representing  the  process  noise. 


Sensor  Measurement 

For  detection  using  an  FM  waveform,  we  only  have  a  target 
positional  measurement  Z^, 


Zk 


%k 

Vk 


+  Vk,  Uk  ~  Af(0,Rk), 


(12) 


where  the  debiased  Cartesian  measurements  (xk,yk)  are 
derived  from  the  time-of-arrival  and  bearing  measurements. 
The  elements  of  measurement  error  covariance  matrix  Rk 
depend  on  the  time  and  bearing  measurement  errors  [9]. 


Kalman  Filter  Equation 

The  target  state  estimate  and  its  covariance  matrix  at  time 
tk  is  denoted  by  X{k\k)  and  P(k\k),  respectively.  Given  the 
initial  target  state  estimate  and  error  covariance  matrix,  the 
recursive  Kalman  Filter  equations  are  described  as  follows: 


X(k\k  -  1)  =  ®k-iX{k  -  1| k  -  1)  (13) 

P{k\k  -  1)  =  $fc_i P(k  -  1| k  -  +  Qk- i  (14) 

L(k)  =  P{k\k  -  1  )CT  ( CP(k\k  -  1  )CT  +  Rk-iY 1  (15) 
X(k\k)=X(k\k-l)+L(k)y(k)  (16) 

P(k\k)  =  {I-L{k)C)P{k\k-l),  (17) 


where  y(k)  =  [xk,  Vk]T  ~  CX(k\k  —  1),  /  is  the  4  x  4  identity 
matrix  and 

r  i  o  o  o 

°  _  0  1  0  0  ■ 

As  for  the  initial  target  state  estimate,  we  use  the  first 
converted  measurement  with  zero  initial  velocity  estimate. 


Tracker  Logic 

The  tracking  algorithm  recursively  processes  the  sets  of 
contacts  from  the  sequence  of  pings  initiated  at  times  (t\  < 
t2  <  •  •  • ).  The  first  ping  at  time  t\  contains  the  set  of  contacts 
{Z3q(l),q  =  I,--  -  ,rij,j  =  I,--  -  Each  of  these  is 

used  to  initiate  a  tentative  track  with  a  state  estimate  and 
its  covariance  matrix.  At  each  subsequent  ping  time  tk ,  we 
have  a  set  of  tracks  (some  of  which  are  confirmed)  from  the 
previous  processing  step  at  time  tk- i,  and  a  current  set  of 
contacts  {ZJq  (k),q  =  1,  •  •  •  ,  rij,  j  =  1,  •  •  •  ,  tir}.  Each  track  is 
predicted  forward  from  time  tk- 1  to  the  current  time  tk .  Then 
the  nearest  neighbor  contact  for  each  track  is  selected  among 
all  validated  contacts.  Validated  contacts  are  those  that  satisfy 
the  following  threshold  condition,  for  a  given  x2  parameter 

2/(fc)T<S'(fc)_1y(fc)  <x2,  (18) 

where  the  innovation  covariance  matrix,  S(k),  is  equal  to 
CP(k\k  —  1  )CT  +  Rk-i,  and  the  nearest  neighbor  is  that 
contact  for  which  the  LHS  of  the  above  equation  is  smallest. 
If  no  validated  contact  exists,  no  track  update  is  performed.  We 
continue  in  this  manner  for  all  tracks,  and  remaining  contacts 


are  used  to  initiate  new  tracks.  A  track  is  confirmed  when  it 
associates  M  contacts  within  N  scans  of  data.  Also,  once  a 
track  is  confirmed,  it  is  terminated  after  K  consecutive  missed 
detections. 


C.  Tracking  Scoring 

The  output  of  the  tracking  algorithm  contains  multiple 
tracks  (tentative  tracks  and  confirmed  tracks).  To  obtain  the 
PT  map,  one  must  map  the  state  estimates  and  their  covariance 
matrices  to  the  probability  of  target  presence  on  the  grid  G. 
Since  not  all  the  tracks  are  originated  from  true  targets,  we 
compute  a  track  score  to  describe  how  likely  (or  unlikely) 
it  is  that  the  reconstructed  track  originates  from  detections 
of  the  true  target  by  considering  kinematic  and  signal-related 
contributions. 

The  track  score  is  defined  as  the  log  of  the  likelihood  ratio 
(LLR)  of  the  hypothesis  that  the  set  of  observations  in  the  track 
are  from  the  same  target  and  the  hypothesis  that  the  track  is  a 
collection  of  false  alarms.  The  initial  track  score  is  estimated 
from  the  a  priori  PT  map  and  the  grid  cell  Pfa.  Thereafter, 
upon  the  receipt  of  data  on  scan  m,  the  score  for  track  n  is 
updated  according  to  the  relationship  [12]: 


Ln{m)  =Ln(m-  1)  +  AL(m),  (19) 


where 
A  L(m) 

A  Lu 

PD 

Pfa 

VC 

S 

d2 


j  ln(l  —  PD)  no  track  update  on  scan  m 
\  A  Ljj  track  update  on  scan  m 


P  D  Vc  }  d 2 

27rPfavl5|  J  2 


estimated  probability  of  detection 
estimated  probability  of  false  alarm 
volume  of  the  measurement  validation  region 
innovation  covariance  matrix 
normalized  statistical  distance  function 


The  track  score  for  track  n  is  updated  recursively  using  Eq. 
(19)  for  all  the  scans  at  each  ping  time.  The  updated  LLR 
score  can  then  be  directly  converted  to  the  probability  of  a 
track  Xj,  being  a  true  target  through  [12]: 


H*i)  = 


exp(Lj) 

1  +  exp (Lj) ' 


(20) 


Let  XA  denote  the  set  of  active  tracks  (tracks  that  are 
updated  using  at  least  one  contact)  at  the  current  ping  time  and 
Na  the  cardinality  of  the  set  XA.  To  convert  the  active  track 
estimates  to  the  probability  on  the  grid  cell,  we  compute  the 
weighting  coefficients  for  those  tracks  in  XA.  When  there  is  no 
contact  updating  a  track  at  the  current  ping  time  (although  the 
track  is  being  propagated  using  KF  propagation  equation),  that 
track  does  not  need  to  be  considered  at  the  current  ping  time 
because  that  track  information  is  not  momentarily  included  in 
evaluation  of  the  positive  PT  map,  but  rather  in  the  negative  PT 
map.  Therefore,  the  track  scores  for  those  tracks  not  included 
in  the  set  XA  are  not  included  in  the  computation  of  the 


844 


weighting  coefficients  or  the  mapping  to  the  positive  PT  map. 
The  weighting  coefficient  for  the  active  track  is  computed  as: 

tin  =  (21) 

D.  Propagation 

A  second  positive  PT  map  is  generated  using  the  KF 
tracker  output.  This  new  positive  PT  map  and  the  negative 
PT  map  from  the  Bayes  update  are  propagated  separately:  the 
positive  map  using  the  KF  tracker  forecasted  estimates  and 
the  negative  map  using  a  diffusion  process  model.  We  assume 
that  the  forecasting  frequency  is  high  relative  to  the  target 
speed  and  that  there  is  no  new  target  at  the  next  forecast  time. 


in  the  x  and  y  directions  are  defined  by  a  normally  distributed 
random  variable  J\f(0,  a2t).  The  Fokker-Planck  (FP)  equation 
for  the  Brownian  motion  is  given  by: 


dp  =  1  2&y  ,  I  2 d2P 
dt  2  dx 2  2  dy 2  ’ 


(25) 


where  p(x,y,t)  represents  the  probability  density  of  target 
presence  at  time  t  in  the  negative  PT  map.  The  FP  equation 
is  solved  using  finite-difference  methods  by  discretizing  the 
variables  x,y,  and  t  as: 


X  - 

-  mS , 

m  =  1,  2,  •  • 

5  x 

(26) 

y  = 

=  nS, 

n  =  1,2,--- 

,  fly 

(27) 

t  - 

=  kAt, 

k  =  1,2,  • 

(28) 

Positive  PT  Map 

The  KF  propagation  equations  (13)  and  (14)  are  used 
to  produce  forecasted  estimates  of  the  active  tracks  at  the 
next  ping  time.  We  then  compute  the  probability  of  target 
presence  on  each  grid  cell  for  each  active  track  using  the 
state  estimate  and  its  covariance  matrix.  Assuming  that  targets 
are  independent  and  the  probability  that  there  is  more  than 
one  target  in  a  grid  cell  is  small,  the  probability  on  each 
grid  cell  accounting  for  all  active  tracks  is  approximated  as 
a  weighted  sum  of  the  individual  probabilities  with  weights 
corresponding  to  their  track  scores.  The  resulting  probability 
map  is  then  normalized  to  the  first  positive  PT  map  that  was 
generated  using  Eq.  (9),  such  that  the  sums  of  all  grid  cell 
probabilities  of  each  map  are  equal.  The  first  positive  PT  map 
is  then  discarded.  We  describe  the  PT  map  conversion  in  the 
following. 

The  probability  of  target  presence  on  each  grid  cell  account¬ 
ing  for  all  active  tracks  is  given  by  a  weighted  sum  of  the 
probability  contributed  by  each  individual  track  as: 

P  (Ti{k  +  1)  =  1|  X(k  +  1|  k),P{k  +  l|fc)) 

NA  , 

p  /  J\f{CXj (k  +  1| k),CPj(k  +  1| k)CT)dxdy, 

7=i  J(x,y)£gi 

(22) 


where  S  is  the  width  of  each  grid  cell  on  G.  The  number  of 
discretized  cells  in  each  direction  is  defined  by  nx  and  ny. 
The  discrete-time  increment  is  defined  by  At  and  the  number 
of  forecast  times  by  nt.  The  above  finite-difference  equations 
are  applied  to  Eq.  (25)  resulting  in: 

a2  o  i 

2  \Pm  —  l,n,k  ^Pm,n,k  + 


Pm,n,k+ 1  —  Pm,n,k  At 


Pm+l,n,fc)  + 


2  Sy2 


(Prn,n—l.k  ^Pm^n^k  Pm,n+l,fc) 


(29) 


Given  the  a  priori  negative  PT  map,  the  above  equation 
updates  the  forecasted  probability  of  target  presence  under  the 
non-detection  hypothesis  as: 


V(Ti(k  +  1)  =  1)  =  Pm,n,k+ 1,  i++(m,n).  (30) 


The  forecasted  PT  map  is  obtained  as  a  summation  of  the 
forecasted  positive  PT  map  and  the  forecasted  negative  PT 
map  given  in  Eqs.  (23)  and  (30)  as: 

P  (Ti(k  +  1)  =  1)  =  P(Tj(fc  +  1)  =  1)+  +  P  (Ti(k  +  1)  =  I)", 

i  =  1,  2,  ■  •  •  ,  Ng. 

(31) 


III.  Simulation 


where  (Xj(k  +  1| k),Pj(k  +  1|&))  denotes  the  state  estimate 
and  covariance  matrix  of  active  track  j,  X(k  +  l\k)  = 
{^(fc  +  llfc),  •  •  •  ,XNA(k  +  l\k)}  and  P(k  +  l\k)  =  {Pi(fc  + 

1| &),■••  ,  Pna  (k+-l\k)}.  The  probability  is  then  renormalized 
to  match  the  sum  of  the  probability  of  the  positive  PT  map 
after  the  Bayes  update  to  obtain  the  forecasted  positive  PT 
map  as: 

P(T;(fc+l)  =  1)+  =  ri-P(Ti(k+l)  =  l\X(k+l\k),P(k+l\k)), 
M  (23) 

iPr(Tt(fc)  =  l|Z(fc))+ 

Efci  P  (Ti(k  +  1)  =  1  \X(k  +  1|A:),  P(k  +  l|fc)) ' 
Negative  PT  Map 

The  negative  PT  map  is  propagated  using  a  Brownian 
motion  diffusion  model  as  described  in  [13].  The  target  motion 


The  probability  of  target  forecaster  is  run  at  the  end  of 
each  ping  to  generate  the  PT  map  at  the  next  ping  time.  In 
this  section,  we  illustrate  that  the  PT  maps  provide  accurate 
multistatic  operational  pictures.  They  include  the  probability 
generated  by  target-originated  contacts  as  well  as  by  false 
contacts  for  further  investigation  in  subsequent  pings.  The 
PT  maps  also  include  the  probability  of  targets  that  might 
be  missed  and  the  probability  that  no  target  is  present. 

In  this  example,  there  are  3  sources  and  4  receivers  deployed 
in  the  search  area  of  40  km  x  40  km.  There  are  two  targets: 
target  1  (Tl)  is  heading  southeast  with  a  constant  speed  of 
6  knots  and  target  2  (T2)  is  heading  east  with  a  constant 
speed  of  5  knots.  The  scenario  is  illustrated  in  Figure  1.  We 
use  a  contact  simulator  to  generate  true  and  false  contacts. 
The  contact  simulator  integrates  a  range-dependent  sonar  per¬ 
formance  model,  an  aspect-dependent  target  strength  model, 
measurement  errors  and  false  contact  distributions  in  number, 
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time-of-arrival,  bearing  and  SNR.  The  simulator  outputs  true 
and  false  contacts.  A  description  of  the  implementation  of  the 
contact  simulator  is  given  in  [10]. 


Field  Layout 


Figure  1.  Simulated  scenario:  three  sources  (red  square),  four  receivers  (green 
circle)  and  two  targets  (blue  line).  Asterisks  indicate  the  start  of  the  tracks. 


dB  to  generate  the  PD  maps.  A  PD  map  for  source  3,  receiver  2 
and  FM  waveform  combination  is  shown  in  Figure  2.  The  track 
initiation  parameters  ( M/N )  are  not  critical  for  this  example 
since  the  forecaster  considers  both  tentative  and  unconfirmed 
tracks.  We  set  the  track  initiation  parameter  (M/N)  to  (1/1) 
and  the  kill  parameter  K  to  8. 


PD  Map,  (S3,R2) 


-20  -15  -10  -5  0  5  10  15  20 


X  (km) 

Figure  2.  The  PD  map  for  source  3,  receiver  2  and  FM  waveform  combination 
generated  with  a  omni-directional  target  strength  of  4  dB  and  a  target  doppler 
of  3  knots. 

We  discretize  the  search  area  into  Ng  =  100  x  100  grid 
cells.  We  begin  the  scenario  with  P(T*(0)  =  1)  =  10-4  for 
i  =  1,2,--*  ,  Ng.  The  Pfa  of  each  grid  cell  is  set  to  0.0012, 
which  corresponds  to  an  average  false  alarm  rate  of  12  per  ping 
interval  in  the  search  area.  In  this  example,  a  ping  is  generated 
every  60  seconds  using  an  FM  waveform  and  a  round-robin 
policy  among  sources.  The  PD  maps  are  precomputed  for  each 
source  and  receiver  pair  using  the  sonar  performance  model. 
We  assume  an  average  omni-directional  target  strength  of  4 
dB,  a  target  doppler  of  3  knots  and  a  signal  fluctuation  of  10 


Figure  3.  The  positive  PT  map  after  Bayes  update  for  ping  1. 
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Figure  4.  The  negative  PT  map  after  Bayes  update  for  ping  1. 

The  positive  and  negative  PT  maps  after  Bayes  update  for 
ping  1  (source  3)  are  shown  in  Figures  3  and  4,  respectively. 
The  PT  values  under  the  detection  hypothesis  increase  in  the 
grid  cells  which  contain  contacts  (Figure  3).  The  PT  values 
under  the  non-detection  hypothesis  decrease  from  the  a  priori 
PT  corresponding  to  the  weighted  average  PD  values  of  the 
grid  cells  (Figure  4). 

Since  many  of  the  contacts  are  false  alarms,  instead  of 
directly  propagating  the  positive  PT  map  using  a  diffusion 
process  model,  the  contacts  are  further  processed  by  the  KF 
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Figure  5.  Track  scores  and  track  weights  for  all  active  tracks  after  ping  1. 
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Figure  6.  Forecasted  PT  map  at  ping  time  2.  White  circles  indicate  the 
locations  of  true  targets. 


tracker  to  obtain  a  more  accurate  representation  of  the  positive 
PT  map.  Figure  5  shows  the  number  of  active  tracks  generated 
after  ping  1  and  their  LLR  scores.  In  this  simulation,  ping  1 
generates  a  total  of  43  contacts  at  the  receivers,  and  there  are 
5  tracks  with  positive  LLR  scores.  Initially,  none  of  the  tracks 
with  positive  LLR  scores  corresponds  to  a  true  target  track. 
The  negative  PT  map  is  propagated  using  a  Brownian  motion 
model  with  a  diffusion  rate  of  5  knots.  The  forecasted  PT 
map  at  ping  time  2  is  obtained  by  combining  the  positive  and 


Figure  7.  Forecasted  PT  map  at  ping  time  2  in  Log  10  scale. 


negative  PT  maps  and  is  shown  in  Figure  6.  The  logarithmic 
scale  graph  for  the  forecasted  PT  map  at  ping  time  2  is  shown 
in  Figure  7  for  better  visualization.  After  ping  1,  PT  values 
in  the  far  east  side  decrease  corresponding  to  the  higher  PD 
values  of  the  grid  cells  in  the  far  east  side.  The  peak  PT  values 
are  distributed  throughout  the  search  area  corresponding  to 
the  contact  locations  with  varying  PT  values  according  to  the 
normalized  track  scores.  Search  effort  must  be  expended  to 
continue  the  true/false  investigation  of  tracks  in  addition  to 
continually  checking  the  grid  cells  that  might  be  overlooked 
due  to  non-detection.  The  forecasted  PT  map  accounts  for  all 
these  probabilities  in  a  single  operational  picture. 

The  number  of  active  tracks  and  their  track  scores  after 
ping  4  are  displayed  in  Figure  8.  After  processing  16  scans 
of  contact  detections,  the  two  active  tracks  which  correspond 
to  the  true  target  tracks  can  clearly  be  seen.  The  forecasted 
PT  maps  at  ping  time  5  in  linear  and  logarithmic  scale  are 
shown  in  Figures  9  and  10.  The  true  target  tracks  correspond 
to  the  highest  PT  values  in  the  forecasted  PT  map.  The  peaks 
corresponding  to  the  false  tracks  are  also  seen  in  the  PT  map. 
The  current  track  weighting  scheme  uses  an  exponential  scale, 
and  therefore  all  the  tracks  with  positive  LLR  are  considered 
as  high-probability  true  tracks  and  are  weighed  nearly  equally. 
The  probability  of  those  tracks  are  continued  to  be  tracked 
under  the  non-detection  hypothesis  even  when  they  are  not 
associated  with  any  contact  in  the  current  ping  (better  seen 
in  Figure  10).  We  can  also  see  that  the  probabilities  in  the 
dark  blue  area  of  Figure  10  decrease  two  order  of  magnitude 
after  4  pings.  This  means  that  the  search  effort  should  now  be 
expended  more  at  the  peaks  and  the  edges.  Different  weighting 
schemes  will  affect  how  much  emphasis  the  search  effort 
places  on  the  active  tracks.  The  tracker  results  can  further 
be  improved  with  a  better  data  association  scheme  such  as 
MHT  [12].  Other  classification  features  can  also  be  included 
for  improved  true/false  discrimination.  All  of  these  aspects 
will  be  addressed  in  future  work. 
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Figure  8.  Track  scores  and  track  weights  for  all  active  tracks  after  ping  4. 
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Figure  10.  Forecasted  PT  map  at  ping  time  5  in  Log  10  scale. 

the  PD  maps  from  a  sonar  performance  model  in  ping  control 
solutions  for  search  missions.  Future  work  will  investigate 
using  the  forecasted  PT  map  as  a  metric  for  ping  control 
optimization.  We  will  also  extend  the  formulation  to  include 
doppler  information  and  explore  different  track  weighting  and 
classification  schemes  for  improved  true/false  discrimination. 
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IV.  Conclusions 

We  have  presented  a  probability  of  target  presence  fore¬ 
caster.  We  have  illustrated  that  the  PT  maps  represent  accurate 
operational  pictures  by  considering  both  detection  and  tracking 
results.  The  forecasted  PT  map  includes  not  only  the  probabili¬ 
ties  of  true/false  tracks  with  emphasis  on  true  tracks  according 
to  their  LLR  track  scores  but  also  the  missed  target  probabil¬ 
ities.  The  forecasted  PT  map  can  be  used  in  conjunction  with 
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